Thermal Expansion

This example is originally from GPUMD and has been added here with minor changes only to demonstrate how to use gpyumd.

1. Introduction

• A given crystal should have a well defined average lattice constant at a given pressure and temperature. Here we use silicon as an example to show how to calculate lattice constants using GPUMD.

Importing Relevant Functions

• The inputs/outputs for GPUMD are processed using the Atomic Simulation Environment (ASE) and the gpyumd package.

import matplotlib.pyplot as plt
import numpy as np
from ase.lattice.cubic import Diamond
from gpyumd.atoms import GpumdAtoms
from gpyumd.io import read_gpumd
import gpyumd.keyword as kwd
from gpyumd.load import load_thermo
from gpyumd.sim import Simulation

2. Preparing the Inputs

• We use a cubic system (of diamond structure) consisting of \(10^3\times 8 = 8000\) silicon atoms and use the minimal Tersoff potential [Fan 2020].

The xyz.in file

Si = GpumdAtoms(Diamond('Si', size=(10,10,10)))
Si.set_max_neighbors(4)
Si.set_cutoff(3)
# We write the xyz.in file here, but gpyumd will also write
# it when we create a simulation.
Si.write_gpumd()
Si

GpumdAtoms(symbols='Si8000', pbc=True, cell=[54.3, 54.3, 54.3])

• The first few lines of the xyz.in file are:

  8000 4 3 0 0 0
  1 1 1 54.3 54.3 54.3
  0 0 0 0 28
  0 0 2.715 2.715 28
  0 2.715 0 2.715 28
  0 2.715 2.715 0 28
  

• Explanations for the first line:

• The first number states that the number of particles is 8000.

• The second number in this line, 4, is good for silicon crystals described by the Tersoff potential because no atom can have more than 4 neighbor atoms in the temperature range studied. Making this number larger only results in more memory usage. If this number is not large enough, GPUMD will give an error message and exit.

• The next number, 3, means the initial cutoff distance for the neighbor list construction is 3 A. Here, we only need to consider the first nearest neighbors. Any number larger than the first nearest neighbor distance and smaller than the second nearest neighbor distance is OK here. Note that we will also not update the neighbor list. There is no such need in this problem.

• The remaining three zeros in the first line mean:

• the box is orthogonal;

• the initial velocities are not contained in this file;

• there is no grouping method defined here.

• Explanations for the second line:

• The first three 1’s mean that all three directions are periodic.

• The remaining three numbers are the box lengths in the three directions. It can be seen that we have used an initial lattice constant of 5.43 A to build the model.

• Starting from the third line, the numbers in the first column are all 0 here, which means that all the atoms are of type 0 (single atom-type system). The next three columns are the initial coordinates of the atoms. The last column gives the masses of the atoms. Here, we show isotopically pure Si-28 crystal, but this Jupyter notebook will generate an xyz.in file using the average of the various isotopes of Si. In some applications, one can consider mass disorder in a flexible way.

Note: We write the xyz.in file here for demonstration purposes, but the Simulation we are about to make also outputs the xyz.in file.

The run.in file

The gpyumd package can be used to generate valid run.in input files as well as other necessary input files. It follows the definitions described in the inputs and outputs documentation for GPUMD.

# Create an empty simulation to build our run.in file with
# Using the correct driver_directory is essential for potential path definitions.
expansion_sim = Simulation(Si, driver_directory='.')

# Programmatically define our thermal expansion runs
pres = 0  # GPa
elast = 53.4059  # GPa
for curr_temp in range(100, 1001, 100):
    curr_run = expansion_sim.add_run(number_of_steps=2e4, run_name=f"T{curr_temp}")

    keywords = [
        kwd.EnsembleNPT.orthogonal(method='npt_ber', initial_temperature=curr_temp,
                               final_temperature=curr_temp, thermostat_coupling=100,
                               barostat_coupling=2000, p_xx=pres, p_yy=pres, p_zz=pres,
                               c_xx=elast, c_yy=elast, c_zz=elast),
        kwd.NeighborOff(),
        kwd.DumpThermo(10)]

    # First run initializations
    if curr_temp == 100:
        keywords.append(kwd.Velocity(initial_temperature=100))
        keywords.append(kwd.TimeStep(dt_in_fs=1))

    for keyword in keywords:
        curr_run.add_keyword(keyword)

# Define the potential
potential_directory = "/path/to/GPUMD/potentials/tersoff"
tersoff_potential = kwd.Potential(filename='Si_Fan_2019.txt', symbols=['Si'], directory=potential_directory)
expansion_sim.add_potential(tersoff_potential)

# Generate xyz.in and run.in files. Can also copy potentials to simulation directory.
expansion_sim.create_simulation(copy_potentials=True)

The run.in input file is given below:

potential Si_Fan_2019.txt 0

velocity 100
time_step 1
ensemble npt_ber 100 100 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 200 200 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 300 300 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 400 400 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 500 500 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 600 600 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 700 700 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 800 800 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 900 900 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

ensemble npt_ber 1000 1000 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 10
run 20000

• The first line uses the potential keyword to define the potential to be used, which is specified in the file Si_Fan_2019.txt.

• The second line uses the velocity keyword and sets the velocities to be initialized with a temperature of 100 K.

• The following 4 lines define the first run. This run will be in the NPT ensemble, using the Berendsen method. The temperature is 100 K and the pressures are zero in all the directions. The coupling constants are 100 and 2000 time steps for the thermostat and the barostat (The elastic constant, or inverse compressibility parameter needed in the barostat is estimated to be 53.4059 GPa; this only needs to be correct up to the order of magnitude.), respectively. The time_step for integration is 1 fs. There are \(2\times 10^4\) steps for this run and the thermodynamic quantities will be output every 10 steps.

• After this run, there are 9 other runs with the same parameters but different target temperatures. Note that the time step only needs to be set once if one wants to use the same time step in the whole simulation. In contrast, one has to use the dump_thermo keyword for each run in order to get outputs for each run. That is, we can say that the time_step keyword is propagating and the dump_thermo keyword is non-propagating.

3. Results and Discussion

• It takes less than 1 min to run this example when a Tesla K40 card is used. The speed of the run is about \(3\times 10^7\) atom x step / second. Using a Tesla P100, the speed is close to \(10^8\) atom x step / second.

Figure Properties

aw = 2
fs = 16
font = {'size'   : fs}
plt.rc('font', **font)
plt.rc('axes' , linewidth=aw)

def set_fig_properties(ax_list):
    ax_list = ax_list if isinstance(ax_list, list) else [ax_list]
    tl = 8
    tw = 2
    tlm = 4

    for ax in ax_list:
        ax.tick_params(which='major', length=tl, width=tw)
        ax.tick_params(which='minor', length=tlm, width=tw)
        ax.tick_params(which='both', axis='both', direction='in', right=True, top=True)

Plot Thermal Expansion

• The output file thermo.out contains many useful data. Here, we load the results and plot the data in the following figure.

thermo = load_thermo()
print("Thermo quantities:", list(thermo.keys()))

Thermo quantities: ['temperature', 'K', 'U', 'Px', 'Py', 'Pz', 'Lx', 'Ly', 'Lz']

time = 0.01*np.arange(1,thermo['temperature'].shape[0]+1) # [ps]
NC = 10  # Number of cells in each direction
NT = 10  # Number of temperature steps
temp = np.arange(100,1001,100)
M = thermo['temperature'].shape[0]//NT
a = (thermo['Lx']+thermo['Ly']+thermo['Lz'])/(3*NC)
Pave = (thermo['Px']+thermo['Py']+thermo['Pz'])/3.
a_ave = a.reshape(NT, M)[:,M//2+1:].mean(axis=1)
fit = np.poly1d(np.polyfit(temp, a_ave, deg=1))

axes = list()
plt.figure(figsize=(12,10))
plt.subplot(2,2,1)
axes.append(plt.gca())
plt.plot(time, thermo['temperature'])
plt.xlim([0, 200])
plt.gca().set_xticks(range(0,201,50))
plt.ylim([0, 1100])
plt.gca().set_yticks(range(0,1101,500))
plt.ylabel('Temperature (K)')
plt.xlabel('Time (ps)')
plt.title('(a)')

plt.subplot(2,2,2)
axes.append(plt.gca())
plt.plot(time, Pave)
plt.xlim([0, 200])
plt.gca().set_xticks(range(0,201,50))
plt.ylim([-0.1, 0.4])
plt.gca().set_yticks(np.arange(-1,5)/10)
plt.ylabel('Pressure (GPa)')
plt.xlabel('Time (ps)')
plt.title('(b)')

plt.subplot(2,2,3)
axes.append(plt.gca())
plt.plot(time, a,linewidth=3)
plt.xlim([0, 200])
plt.gca().set_xticks(range(0,201,50))
plt.ylim([5.43, 5.48])
plt.gca().set_yticks([5.44,5.46,5.48])
plt.ylabel(r'a ($\AA$)')
plt.xlabel('Time (ps)')
plt.title('(c)')

plt.subplot(2,2,4)
axes.append(plt.gca())
Tpoly = [0, 1100]
plt.plot(Tpoly, fit(Tpoly),color='C3')
plt.scatter(temp, a_ave,s=200,zorder=100,facecolor='none',edgecolors='C0',linewidths=3)
plt.xlim([0, 1100])
plt.gca().set_xticks(range(0,1101,500))
plt.ylim([5.43, 5.48])
plt.gca().set_yticks([5.44,5.46,5.48])
plt.ylabel(r'a ($\AA$)')
plt.xlabel('Temperature (K)')
plt.title('(d)')

set_fig_properties(axes)
plt.tight_layout()
plt.show()

(a) Instant temperature as a function of simulation time. (b) Instant pressure as a function of simulation time. (c) Instant lattice constant as a function of simulation time. (d) Averaged lattice constant (over the last 10 ps in each run for a given temperature) as a function of temperature.

Additional figure notes: - (a): The temperature for each run quickly reaches the target temperature (with fluctuations).

• (b): The pressure (averaged over the three directions) for each run quickly reaches the target pressure zero (with fluctuations).

• (c): The lattice constant (averaged over the three directions) for each run reaches a plateau (with fluctuations) after some steps.

• (d): We calculate the average lattice constant at each temperature by averaging the second half of the data for each run. The average lattice constants at different temperatures can be well fit by a linear function, with the thermal expansion coefficient being estimated to be $:nbsphinx-math:alpha `:nbsphinx-math:approx 7.2 :nbsphinx-math:times `10^{-6} $ K-1.

4. References

[Fan 2020] Zheyong Fan, Yanzhou Wang, Xiaokun Gu, Ping Qian, Yanjing Su, and Tapio Ala-Nissila, A minimal Tersoff potential for diamond silicon with improved descriptions of elastic and phonon transport properties, J. Phys.: Condens. Matter 32 135901 (2020).